function plotrocks(filenumber)
S=readdata(0,2); T=readdata(0,3);
dx=readdata(0,8);
dz=readdata(0,9);
x=zeros(1,S); z=zeros(1,T);
for s=1:S-1 x(s+1)=x(s)+dx(s); end
for t=1:T-1 z(t+1)=z(t)+dz(t); end
depthrange=z;

temphist=zeros(S,length(filenumber),length(depthrange));
Phist=zeros(S,length(filenumber),length(depthrange));
refsurface=readdata(0,25); refsurface=refsurface(:,4);
surfacehist=zeros(S,length(filenumber)); avtemphist=surfacehist; timehist=zeros(1,length(filenumber));
timeindex=1;

for count=filenumber
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
    %     multigridlevels=readdata(count,4);
    %     bound=readdata(count,5);
    time=readdata(count,6);
    %     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
    %     mcs=readdata(count,12);
    %     mct=readdata(count,13);
    %    mtype=readdata(count,14);
    mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
    %     vx=readdata(count,21);
    %     vz=readdata(count,22);
    P=readdata(count,23);
    Nsurfaces=readdata(count,24);
    surfaces=readdata(count,25);

    count
    x=zeros(1,S); z=zeros(1,T);
    for s=1:S-1 x(s+1)=x(s)+dx(s); end
    for t=1:T-1 z(t+1)=z(t)+dz(t); end
    xp=zeros(1,S-1); zp=zeros(1,T-1);
    for s=1:S-1
        xp(s)=0.5*(x(s)+x(s+1));
    end
    for t=1:T-1
        zp(t)=0.5*(z(t)+z(t+1));
    end
    
    [Zp,Xp] = meshgrid(zp,xp); P=P(1:end-1,1:end-1);
    [gridtemp,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1);
    [gridP,xgrid,zgrid]=bilingrid(P(:),Xp(:),Zp(:),dx,dz,1,1);
   
    [Zgrid,Xgrid] = meshgrid(zgrid,xgrid);
    for d=1:length(depthrange)
        %localtemp=interp2(Xgrid',Zgrid',gridtemp',x,depthrange(d));
        %localtemp=gridtemp(:,d);
        temphist(:,timeindex,d)=gridtemp(:,d);
        Phist(:,timeindex,d)=gridP(:,d);
    end
    integrationdepths=linspace(40e3+0e3,max(z),2*T);
%     for s=1:S
%         localtemp=interp2(Xgrid',Zgrid',gridtemp',x(s),integrationdepths);
%         avtemphist(s,timeindex)=trapz(integrationdepths,localtemp);
%     end
    
    sealevel=0.5*(surfaces(1,Nsurfaces)+surfaces(S,Nsurfaces));
    %surfacehist(:,timeindex)=surfaces(:,4);
    %surfacehist(:,timeindex)=surfaces(:,4)-refsurface;
    tempsurface=surfaces(:,4)-0*refsurface; surfacehist(:,timeindex)=tempsurface-mean(tempsurface);  %correction for long term drift
    timehist(timeindex)=time/(60*60*24*365*1e6);
    timeindex=timeindex+1;
end
save surface_time timehist surfacehist temphist avtemphist depthrange Phist x;